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Abstract 

This paper revisits the fundamental equations for the solution of the frictionless uni¬ 
lateral normal contact problem between a rough rigid surface and a linear elastic 
half-plane using the boundary element method (BEM). After recasting the result¬ 
ing Linear Complementarity Problem (LCP) as a convex quadratic program (QP) 
with nonnegative constraints, different optimization algorithms are compared for 
its solution: (i) a Greedy method, based on different solvers for the unconstrained 
linear system (Conjugate Gradient CG, Gauss-Seidel, Cholesky factorization), (m) 
a constrained CG algorithm, (Hi) the Alternating Direction Method of Multipliers 
(ADMM), and (iv) the Non-Negative Least Squares (NNLS) algorithm, possibly 
warm-started by accelerated gradient projection steps or taking advantage of a load¬ 
ing history. The latter method is two orders of magnitude faster than the Greedy 
CG method and one order of magnitude faster than the constrained GG algorithm. 
Finally, we propose another type of warm start based on a refined criterion for the 
identification of the initial trial contact domain that can be used in conjunction 
with all the previous optimization algorithms. This method, called Gascade Multi- 
Resolution (CMR), takes advantage of physical considerations regarding the scaling 
of the contact predictions by changing the surface resolution. The method is very 
efficient and accurate when applied to real or numerically generated rough surfaces, 
provided that their power spectral density function is of power-law type, as in case 
of self-similar fractal surfaces. 
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1 Introduction 


Contact mechanics between rough surfaces is a very active area of research for 
its paramount importance to address several practical applications in physics 
and engineering. Understanding the evolution of the contact domain and con¬ 
tact variables, such as load, real contact area, contact stiffness, and many 
others, that depend on the morphological properties of roughness, is still con¬ 
sidered a challenging problem today. The reader is referred to (Barber, 2003; 
Nosonovsky and Bhushan, 2005; Ciavarella et ah, 2006; Hyun and Robbins, 

2007; Ciavarella et ah, 2008a,b; Carbone and Bottiglione, 2008; Paggi and Ciavarella, 
2010; Campaha et ah, 2011; Paggi and Barber, 2011; Paggi et ah, 2014; Yastrebov, 
2014) for an overview of research results developed during the last decade. 

Semi-analytical contact theories that are able to provide synthetic predic¬ 
tions of the contact response is also a challenging topic. A comparison and 
validation on benchmark results is necessary to understand the limitations 
of existing approaches and propose further improvements. Experimental in¬ 
vestigations are difficult to make and involve approximations, for example 
very often the contact parameters can only be estimated by indirect measure¬ 
ments of thermal or electric resistances of compressed rough joints (McCool, 

1986; Sridhar and Yovanovich, 1994) or are mostly limited to measurements 
of real contact area under special conditions (O’Callaghan and Probert, 2005; 
Hendriks and Visscher, 1995). Therefore, numerical methods are essential to 
acquire as much information as possible about the contact problem at hand 
and infer general conclusions. 

In spite of its effectiveness and versatility, the hnite element method (FEM) 
has been mainly applied in mechanics to solve contact problems between 
rough surfaces in which the constitutive behavior of the bulk is not linear 
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elastic. For instance, the study of elasto-plastic contact problems with rough¬ 
ness (Hyun et ah, 2004), where an explicit approach was used to reduce the 
high computational cost, and the study involving frictional dissipative phe¬ 
nomena in visco-elastic materials, where the energy dissipation in the bulk is 
essential and can be well predicted by FEM (Wriggers and Reinelt, 2009), are 
worth mentioning. 

In the linear elastic regime, when the multi-scale character of roughness cover¬ 
ing a wide spectrum of wavelengths is the main focus, the use of the boundary 
element method (BEM) is historically preferred over FEM (Andersson, 1981; 
Man, 1994). This is essentially due to the fact that only the surface must be 
discretized and not the bulk. Moreover, it is not necessary to adopt surface 
interpolation techniques, like Bezier curves, to discretize the interface (see, 
e.g., the rigorous studies in (Wriggers, 2006, Ch. 9) and (Hyun et ah, 2004)), 
which must be used with care to avoid smoothing out artificially the fine scale 
geometrical features of roughness. 

In the application of BEM, the frictionless contact problem between two lin¬ 
ear elastic rough surfaces is mathematically equivalent to the problem of the 
normal contact between a rigid rough surface and an elastic half-plane with 
equivalent elastic parameters, see (Barber, 2003) for a rigorous proof. The 
core of BEM is based on the so-called Green’s functions, that relate the dis¬ 
placement of a generic point of the half-plane to the action of a concentrated 
force on the surface caused by contact interactions. An integral convolution 
of all the contact tractions provides the deformed contact conhguration. After 
introducing a discretization of the half-plane consisting of a grid of boundary 
elements, the problem of point-force singularity is solved numerically by us¬ 
ing the closed-form solution for a patch load acting on a finite-size boundary 
element (Johnson, 1985, Ch. 3,4). The contact problem is then set in terms of 
equalities and inequalities stemming from the unilateral contact constraints 
and can be solved by constrained optimization. In this regard, apart from the 
discretization error intrinsic in any numerical method, BEM provides the high¬ 
est attainable accuracy for discrete problems (Polonsky and Keer, 1999). The 
basic version of BEM can be also extended to solve rough contact problems 
with friction (Li and Berger, 2003; Pohrt and Li, 2014) and between viscoelas- 
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tic materials (Carbone and Putignano, 2013). 


With the aim of investigating the effect of roughness at multiple scales, the 
availability of computational methods that can solve large contact problems 
in an efficient and fast way is of crucial importance. The size of the linear 
system of equations relating the contact pressures to the normal deflections 
can be in fact quite large, as it arises from high resolution profilometric surface 
samples of 512x512 heights and very large indentations. Hence, the compu¬ 
tational challenges regard two main aspects: (i) efficiently solve the system 
of linear equations; {ii) impose the satisfaction of the unilateral contact con¬ 
straints (contact inequalities). Regarding the first issue, iterative methods like 
the Conjugate Gradient algorithm or the Gauss-Seidel method (Francis, 1983; 
Borri-Brunetto et ah, 1999, 2001) have been widely used. Alternatively, the 
capabilities of multigrid or multilevel methods have been exploited (Raous, 
1999; Polonsky and Keer, 1999) to approximately solve the equation system on 
coarse grids and then project the results on finer grids. Finally, we mention the 
fast method and its variants based on the solution of the linear system of equa¬ 
tions in the Fourier space (see, e.g., (Nogi and Kato, 1997; Polonsky and Keer, 
2000a,b; Batrouni et ah, 2002; Scaraggi et ah, 2013; Prodanov et ah, 2014)). 

Regarding the imposition of the contact inequalities, (Johnson, 1985, p.l49- 
150) suggested to apply a greedy approach: after solving the equation set for 
the unknown tractions, the boundary elements for which these are negative 
(tensile) are excluded in a following iteration from the assumed contact area 
and the corresponding pressures set equal to zero. Johnson (1985) [p.149-150] 
stated that “experience conhrms that repeated iterations converge to a set of 
values of pressures which are positive where contact takes place and zero oth¬ 
erwise”. To the best of the authors’ knowledge, a rigorous proof of convergence 
of this method has not been found in the literature. However, if valid, it allows 
to use any numerical method to solve the unconstrained set of linear equa¬ 
tions and then impose a correction in a recursive way. Indeed, this numerical 
approach has been successfully applied by many authors, such as Kubo et ah 
(1981) and Borri-Brunetto et ah (1999, 2001) who used this greedy approach 
in conjunction with a Gauss-Seidel iterative algorithm for the solution of the 
unconstrained set of linear equations, and Karpenko and Akay (2001) and 
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Batrouni et al. (2002) who applied it together with a numerical scheme based 
on the Fast Fourier Transform (FFT). 

In alternative to the greedy approach, Polonsky and Keer (1999) proposed 
a constrained Conjugate Gradient method based on the theory in (Hestenes, 
1980, Ch. 2,3) to solve the linear system of equations and rigorously impose the 
satisfaction of the contact constraints. For the solution of the system of equa¬ 
tions, a multi-grid solution scheme was proposed in (Polonsky and Keer, 1999) 
and then a FFT algorithm was considered in (Polonsky and Keer, 2000a,b). 

In this paper, we first examine the validity of the greedy approach based on 
a monotonic elimination of tensile points. We show that this approach usu¬ 
ally hnds the exact solution but, as we prove by a counter-example, it may 
fail. Then, we show that other optimization algorithms such as Non-Negative 
Least Squares (NNLS) and the Alternative Direction Method of Multipliers 
(ADMM) can be used in alternative to the greedy approach, by exploiting the 
equivalence between the contact problem and quadratic programming with 
unilateral non-negativity constraints. Moreover, we propose warm starting 
techniques for the optimization algorithms that are especially useful in case 
of a solution of a sequence of increasing or decreasing displacements. 

This paper provides a comprehensive comparison of the computational per¬ 
formance of the greedy approach (used in conjunction with different uncon¬ 
strained solvers like the Conjugate Gradient, the Gauss-Seidel iterative scheme, 
or the MATLAB’s mldivide solver ^), of the original constrained CG method 
by (Polonsky and Keer, 1999), and of novel optimization algorithms that are 
able to exploit warm starts for solving convex quadratic programs subject 
to non-negativity constraints. As a main conclusion, the proposed NNLS al¬ 
gorithm with warm start based on accelerated gradient projections (GPs) is 
found to be one order of magnitude faster than the algorithm by Polonsky and Keer 
(1999) and two orders of magnitude faster than the greedy approach. 

Finally, by exploiting the morphological features of the contact domain of 

^ According to documentation, mldivide solves linear systems with sym¬ 
metric positive definite matrices by computing a Cholesky factorization, see 
http://www.mathworks.it/help/MATLAB/ref/mldivide.html 
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fractal surfaces, we propose in this paper a cascade multi-resolution algorithm 
that can further reduce computation time by at least a factor two with respect 
to the NNLS algorithm with accelerated GPs. 


2 Mathematical formulation 

In the framework of BEM, the normal displacements m(x) at any point of 
the half-plane identihed by the position vector x are related to the contact 
pressures p{y) at other points as follows (Johnson, 1985; Barber, 2010): 



( 1 ) 


where iJ(x, y) represents the displacement at a point x due to a surface con¬ 
tact pressure p acting at y and S is the elastic half-plane. For homogeneous, 
isotropic, linear elastic materials, the influence coefficients are: 



( 2 ) 


where E and u denote, respectively, the composite Young’s modulus and Pois¬ 
son’s ratio of half-plane, and || • || the standard Euclidean norm. The total 
contact force P is the integral of the contact pressure held 



( 3 ) 


By referring to Fig. 1, in the following we dehne for each surface point x G S' 
its elevation ^(x), measured with respect to a reference frame, and set ^max — 
maXxgs ^(x) the maximum elevation. The indentation of the half plane at the 
points in contact is denoted by it, whereas a generic displacement along the 
surface is u. 

We consider the following problem: 

Problem 1 For a given far-field displacement A > 0 in the direction perpen¬ 
dicular to the undeformed half-plane, find the solution of the normal contact 
problem n(x), p(x) satisfying (1) and the unilateral contact (linear comple- 
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(C) 



Figure 1. Sketch of the contact problem between a rigid rough surface and an elas¬ 
tic half-plane. Its deformed configuration corresponding to the imposed far-held 
displacement A is depicted with a black solid line. The red dashed line correspond¬ 
ing to a rigid-body motion of the half-plane identihes the heights to be included in 
the initial trial contact domain. Once Problem 1 is solved we may have: (i) heights 
certainly not in contact from the beginning, type (a); (ii) heights loosing contact 
due to elastic interactions, type (b); {iii) heights in contact, type (c). 

mentarity) conditions 

m(x) — ■u(x. A) > 0, (4a) 

p(x) > 0, (4b) 

(m(x) - -(/(x, A))p(x) = 0, (4c) 

for all points x G S', where contact tractions are positive when compressive. 

Introducing the quantity ta(x. A) = m(x) — -u(x. A), Eq.(4) can be rewritten 
as: 


tc(x. A) > 0, 

(5a) 

p(x) > 0, 

(5b) 

tc(x, A)p(x) = 0. 

(5c) 


Problem 1 is an infinite-dimensional linear complementarity problem. We find 
a finite-dimensional approximate solution by discretizing the surface as a 
square grid of spacing 5 consisting oi N xN average heights. Let S'^ be the cell 
of area 5^ indexed by i, j G /at, with 1^ — {1, x {1,..., A}. Let Xjj = 

L&Sij xdx, = /xe 5 ,, '^(x)dx, pij A 4g^,^p(x)dx, and A M(x)dx 
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be, respectively, the barycentric coordinate, average height, resultant of the 
contact pressures, and the corresponding displacement on the surface element 
Sij. Consider the following discretized version of (1) 

N N 

~ ^i-k,j-lPk,l ( 6 ) 

k=l 1=1 

for all {i,j) G In, where the term is the Green function used in (1) 

averaged over the elementary area 5^, which corresponds to the displacement 
induced by a uniformly loaded square: 

= [ if(x,y)dydx, (7) 

0 J Sij ^ S]zi 

and 


Vk,i > 0, V(fc,/) e In- 


( 8 ) 


For instance, Borri-Brunetto et ah (1999) used the following approximation 
related to a uniform pressure acting on a rounded patch of radius 5/2: 


Hi—kA — l ' 


EttS^ 

2 


arcsm 


2||xjj - Xfc,i 


ii i = k and j = I 
a i ^ k, j ^ I 


( 9 ) 


but other formulae for a square patch can also be taken as in (Pohrt and Li, 
2014). 


Let Ic = {{hi) ^ In ■ ^i,j < Cmax — be the set of indices corresponding to 
elements Sij that are certainly not in contact (cf. Fig. 1), and hence 


Pk,i — 0,V(/i:,/) G Ic, (10) 

let m = 4^^Ic be the number of elements of Ic and n = ^Ic the number of 
elements belonging to the initial trial contact domain, Ic = In\Ic- The set Ic 
is only a superset of the set //. of actual contact points, since the corrections to 
the displacements induced by elastic interactions may induce lack of contact 
in some elements {i,j), i.e., Uij > Uij, where hjj = A — ^max + 6j is the value 
of the compenetration of the height corresponding to the element {i,j) in the 
half-plane (see Fig. 1). 







For a generic {i,j) G Ic corresponding to an element of the surface which is 
potentially in contact with the elastic half-plane, we denote by 

^0 ( 11 ) 

the corresponding elastic correction to the displacement. Clearly, it must hold 
that 


= 0, V(h j) e Ic (12) 

since Wi^ > 0 implies no contact between the surfaces and therefore no pres¬ 
sure, while Pij > 0 implies contact, Uij = Uij, or equivalently Wij = 0. 

By taking into account that pk,i = 0 for all {k, 1) G Iq, Eq. (6) can be recast 
as the following condition 

Wi,j + Ui,j= X) Hi_kj-iPk,i, e Ic, (13) 

{k,lWc 

which is limited to the nodes belonging to the initial trial contact domain Ic, 
whose number of elements is in general signihcantly smaller than those of In- 
The relations (8)-(13) can be recast in matrix form as the following Linear 
Complementarity Problem (LCP) (Cottle et ah, 1992): 

w = Hp — u (14a) 

w > 0, p > 0, w'p = 0, (14b) 

where w G is the vector of unknown elastic corrections Wij, {i,j) e Ic, 
w' denotes its transpose, p G M" is the vector of unknown average con¬ 
tact forces pij, {i,j) G Ic, u G M" is the vector of compenetrations 
{i,i) G Ic, and H = H' is the matrix obtained by collecting the compli¬ 
ance coefficients Hi_kj-i, for {i,j), {k,l) G Ic- Due to the properties of linear 
elasticity (Johnson, 1985, p.l44) we have that 

H = H' ^ 0, (15) 

that is H is a symmetric positive dehnite matrix (with the additional property 
deriving from (9) of having all its entries positive). After solving (14), the 
vector u G M” of normal displacements Uij, {i,j) G Ic, is then simply retrieved 
as u = u -|- w. 
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By the positive definiteness property (15) of H, we inherit immediately the 
following important property (Cottle et ah, 1992, Th. 3.3.7): 

Property 1 The discretized version (8), (10)-(13) of Problem 1 admits a 
unique solution p, u, for all A > 0. 


The TCP problem (14) corresponds to the Karush-Kuhn-Tucker (KKT) con¬ 
ditions for optimality of the following convex quadratic program (QP) 


miup 2P “ ^'P 

s.t. p > 0 


(16a) 

(16b) 


in that the solution p of (16) and its corresponding optimal dual solution w 
solve (14), and vice versa. 


Problem (16) is consistent with former pioneering considerations by Kalker and van Randen 
(1972) and also summarized in (Johnson, 1985, p. 151-152). In fact, the contact 
pressures solving the unilateral contact problem can be obtained by minimiz¬ 
ing the total complementary energy W of the linear elastic system, subject to 
the constraint p(x) > 0, Vx G S. For a continuous system, the total comple¬ 
mentary energy is 


W 


U — [ p(x)m(x. A) dx, 
Js 


(17) 


where U is the internal complementary energy of the deformed half-plane in 
contact. For linear elastic materials, we have: 

U = ^^p(x)M(x)dx. 



Although such an energy-based approach can be used to derive FEM for¬ 
mulations, it is also possible to remain within BEM and introduce a surface 
discretization as before. By invoking the averaged Green’s functions in (7), 
the discretized version of W leads to 


^ o ^ ^ ^ pijUij (19) 

hj)6/c {k,l)&Ic {iJWc 

which represents a quadratic function of p to be minimized, under the con¬ 
straints pij > 0, V(i, j) G Jc, as in (16). Since it is unlikely that the contact 
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area is known a priori, the active set of nodes in contact results only after 
solving problem (14) or equivalently (16). 

A large variety of solvers for LCP and QP problems were developed in the last 
60 years (Beale, 1955; Fletcher, 1971; Goldfarb and Idnani, 1983; Cottle et ah, 
1992; Schmid and Biegler, 1994; Patrinos and Bemporad, 2014), and is still 
an active area of research in the optimization and control communities. His¬ 
torically, in the mechanics community, Kalker and van Randen (1972) pro¬ 
posed the simplex method, although it was found to be practical only for 
relatively small N. More recent contributions adopt algorithms to solve the 
unconstrained linear system of equations and then correct the solution by 
eliminating the boundary elements bearing tensile tractions (Francis, 1983; 
Borri-Brunetto et al., 1999, 2001), or use a constrained version of the Con¬ 
jugate Gradient (CG) algorithm (Polonsky and Keer, 1999). These methods 
are simply initialized by considering arbitrary nonnegative entries in p, with¬ 
out taking advantage of the monotonic increase (or decrease) of pressures by 
increasing (or decreasing) the far-held displacement, an important property 
guaranteed by rigorous elasticity theorems (Barber, 1974). The history of pres¬ 
sures is saved during a contact simulation and it is easy to access and use and 
it can be benehcial to save computation time. 

Next section presents effective optimization algorithms for solving the QP 
problem (16) and compares their performance with respect to the Greedy 
CG method. Contrary to the latter, not only the considered QP have the 
guaranteed property of always converging to the unique solution p, u for any 
given A > 0, but also the history of loading can be more efficiently taken into 
account as a warm-start, with a signihcant saving of computation time. 


3 Optimization algorithms 


Since now on, we use the subscript i to denote the i-th component of a vector 
or the i-th row of a matrix, the subscript X to denote the subvector obtained 
by collecting all the components i G X of a vector (or all the rows i of a 
matrix), and the double subscript X,Xi to denote the submatrix obtained by 
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collecting the i-th row and j-th column, for all i G X, j G Xi. 


3.1 Greedy methods 


A greedy method corresponds to solve problem (16) by iteratively solving the 
unconstrained linear system of equations w = Hp — u = 0 with respect to 
p and increasingly zeroing negative elements of p until the condition p > 0 
is satished. By construction we obtain w'p = 0. The method is described in 
Algorithm 1, in which a standard Conjugate Gradient employed to solve the 
unconstrained linear system of equations. Steps 2.1-2.4 can be replaced by 
any other algorithm for solving the linear system of equations, like the Gauss- 
Seidel iterative scheme as in (Borri-Brunetto et ah, 1999, 2001), the MAT- 
LAB’s mldivide solver, or even the FFT algorithm as in (Karpenko and Akay, 
2001; Batrouni et ah, 2002). 

Assuming that the prescribed initial p and X are such that pj = 0 for all 
i G {1,... ,n} \ X, and iF^ax is sufficiently large, the output of the greedy 
algorithm leads to a contact pressure vector p* and a normal displacement 
vector u* satisfying u* = Hp*, p* > 0, (u* — u)'p* = 0. In fact, condition 
p* > 0 is guaranteed by the condition in Step 2 up to e precision. By letting 
w* = u* — u, at termination of the algorithm we have Wj = Hj ,xpj - ui = 0 
because of the solution of the CG method (Step 2.4), or equivalently Uj = Uj 
(cf. Step 5). By setting uj = H f jpx in Step 5, and recalling that pJ = 0, we 
have 


Wi 


0 0 


Px 


0 


Hi,i Hx x 


Px 


-Ux 






+ 






+ 




X 0 


0 




Hx,x Hx^x 




-Ux 


and hence u* = w* -|- u = Hp*. The complementarity condition (u* — u)'p* = 
(w*)'p* = 0 follows by construction, as Step 2.4 zeroes all the components of 
w* that correspond to nonnegative p*, Vj G X, and zeroes all the components 
p* that correspond to possible nonzero components w*, Vj G X. 

However, to the best of the authors’ knowledge, no formal proof exists that the 
condition w|- > 0 is satished after the algorithm terminates, i.e., that u* > u. 


12 
















Algorithm 1. Greedy method with Conjugate Gradient (greedy CG) 


Input: Matrix H = H' >- 0, vector 
X C {1,..., n} such that P{i,,..,n}\i = 

tolerance e > 0. 

u; initial guess p and initial active set 

0; maximum number iCmax of iterations. 

(1) f ^ 0; X 

t— {1,..., n} \ X; 



(2) while {i 

< JPmax and min(p) < 

— e) or f = 

0 do: 

(2.1) wx Hj^xpx — Uj; 



(2.2) Uu, ||wx||2; 



(2.3) bj < - Wx 



(2.4) while > e and i < iCmax 

do: 


(2.4.1) 

sx ^ Hx,xbx; 



(2.4.2) 

PI ^ PI bi; 



(2.4.3) 

Wx t— Hx,xPi ~ hx; 



(2.4.4) 

bx^ wx+bjsxbi; 



(2.4.5) 

Wx ^ Wx; 



(2.4.6) 

riyo ^ llwxlh; 



(2.4.7) 

z •<— i + 1; 



(2.5) for j EX do: 



(2.6.1) 

if pj < —e then pj i 

0; X^X\ 

{j}-, i ^ ^ {j}; 

(3) p* ^ p; 




(4) U*j- = Ux: 

, uj Hx^xPi; 



(5) end. 




Output: Contact force vector p* and normal displacement vector u*. 


If the algorithm is applied to randomly generated u vectors and H positive 
dehnite matrices with positive coefficients, in many cases the LCP is not solved 
exactly. In contact mechanics, the only evidence that this condition is satished 
has been shown in simnlations (see, e.g., (Batronni et ah, 2002)). Indeed, we 
obtained the following connterexample in which the greedy method failed in 
getting the solntion also for H whose coefficients are given by Eq.(9) ^. 

Example 1 Consider a sqnare mesh with grid spacing 6 consisting of iV x iV 
bonndary elements indexed by (i, j) G In, In = , N} x {1, • • •, N}- Snp- 

^ The MATLAB routine of the counterexample is available for download at 
http://musam.imtlucca.it/counterexample.m 
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pose that all the boundary elements are included in the initial trial contact do¬ 
main Ic {n = N X N) by assigning a positive value of Ujj to all elements. This 
may represent a situation where a cluster of densely packed heights comes into 
contact. Since Uij depends on the height held which is a random variable, 
for the sake of generality we extract the values Ujj- randomly from a uniform 
distribution in the interval (0,1). The matrix H is assembled according to (9). 
By running a sequence of 100 random simulations, we hnd that in approxi¬ 
mately 40% of the simulations the greedy method provides a solution which 
violates the condition w*j > 0 in at least one element. This lack of getting 
the right solution is observed for any size n of the problem. One of the wrong 
results obtained for n = 100 is shown in Fig. 2. The assigned random values 
of u are plotted in Fig. 2(a) for the sequence of boundary elements (from 1 
to 100) composing the mesh. The solution w* presents a negative entry in 
one single element (element 62 in Fig. 2(b)). The computed contact forces are 
compared in Fig. 2(c) with the values corresponding to the exact solution of 
the problem (green dots) obtained by using the NNLS algorithm presented 
in Section 3.3, that is proven to satisfy the LCP conditions (14) exactly. Al¬ 
though just one value of w* is negative, the overall solution is affected by this 
violation. We observe in fact a false contact detection for the element number 
62 violating the condition w*^- > 0, a contact not detected (element 81) and 7 
contact forces signihcantly underestimated with respect to the exact ones. ■ 


For less densely packed boundary elements belonging to Ic. for instance with 
a minimum distance of 25 between them instead of 5 as in Example 1, the 
algorithm was found to always provide a solution satisfying the condition 
w* > 0. Other benchmark tests considering a deterministic smooth variation 
of u, as in case of an indentation by a smooth sphere or by a flat punch, did not 
show any convergence problem to the solution as well, although the boundary 
elements in contact are densely packed as in the counterexample shown before. 
In conclusion, although it is likely that the diagonally dominant property of 
the matrix H plays a role in the robustness of the algorithm, it remains an open 
problem to hnd exact mathematical requirements for H and u that guarantee 
the greedy method to provide a solution satisfying w* > 0, so that all the 
LCP conditions (14) are met. 
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(a) Random u as input 



(b) Computed wj 



Boundary element 
(c) Computed p/E 

Figure 2. Counterexample showing that the greedy CG method fails in getting 
the correct solution (5 = 1 a.u. of L, as u and w*; E = 0.01 F/L"^). Green dots 
correspond to the correct contact forces satisfying the LCP and are obtained by 
using the NNLS method, Sec. 3.2 
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Therefore, as a word of caution, the reliability of the greedy method should be 
carefully checked in case of applications of BEM to contact problems governed 
by other forms of H, as in the case of contact with an anisotropic or an 
inhomogeneous half-space, or in the presence of multiple helds. 

Another drawback of the algorithm is the difficulty to warm start the method 
with a proper choice of the initial active set X. Since at Step 2.5.1 the number 
of elements in the sequence X is decreased by removing negative enough com¬ 
ponents Pj of the current solution vector, i.e., eliminating the points bearing 
tensile (negative) forces, in a monotonic way (no index j that has been removed 
from X can be added back), a safe cold start is to set X = {1,..., n} and pick 
up a vector p > 0, usually a vector with arbitrary non-negative numbers. The 
history of contact forces obtained during the solution of a sequence of imposed 
displacements is not taken into account by the method to accelerate its con¬ 
vergence, although we know that contact forces are monotonically increasing 
functions of the far-held displacement. In any case, for a complex sequence of 
loading with an increased or decreased far-held displacement, any warm start¬ 
ing on forces cannot be implemented in the method, since the elimination of 
contact points is irreversible. 


3.2 Constrained Conjugate Cradient 


A constrained CG algorithm was proposed by Polonsky and Keer (1999) based 
on the theory by (Hestenes, 1980, Ch. 2,3) to solve the linear system of equa¬ 
tions and rigorously impose the satisfaction of the contact constraints. Algo¬ 
rithm 2 has been applied by Polonsky and Keer (1999) to simulations under 
load control. However, it can be used also for displacement control. The con¬ 
dition for convergence set by Polonsky and Keer (1999) in terms of relative 
variation in the local contact forces from an iteration to the next has been re¬ 
cast in terms of the error in the local contact displacements. The two criteria 
are completely equivalent. 

This constrained CG algorithm does not remove the points bearing tensile 
forces from the active set. Therefore, the size of the linear system of equations 
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Algorithm 2. Constrained Conjugate Gradient 


Input: Matrix H = H' >- 0, vector u, initial guess p > 0, initial active set 
X = {1,..., n}; maximum number iCmax of iterations, tolerance e > 0. 


(1) i ^ 0, n^,oid = 1, ci = 0, err = +cx); 

(2) w ^ Hp - ii; 

(3) while {i < /fmax and err > e): 

Tl 

(3.1) if i = 0 then t ^ w else: t w + d- 


(3.2) r = 

V ^ +/1 


^td,old 


^old) 


t'Ht’ 

(3.3) p ^ p — rt; 

(3.4) Vj G X : pj G- max{pj, 0}; 

(3.5) Find loi = {j G X : pj = 0, Wj < 0}; 

if IqI 0 then d 1 else d 0, pj ■i pj TXVj^ Vy G 

(3.6) X G- {j : Pj > 0} U lau 

(3.7) fold ^ fj ^Jiijold '^wi 

(3.8) w ^ Hp — u; 

(3.9) = ||w||2; 

(3.10) err ^ \n^ - n^,oid|/?^«;,oid; 

(3.11) + 

(4) p* ^ p; u* = Hp*; 

(5) end. 

Output: Contact force vector p* and normal displacement vector u*. 


is not reduced during the iterations, increasing the computation time for its 
solution. On the other hand, the method assures the satisfaction of the LCP 
conditions (14) and it is found to convergence with a reduced number of iter¬ 
ations as compared to the Greedy CG algorithm. Although not investigated 
in (Polonsky and Keer, 1999), it can be warm started in case of a sequence 
of loading steps by considering both an initial trial contact domain and a set 
of contact pressures derived from the previous converged solution. The FFT 
method can be used to accelerate step (2.7) as in (Polonsky and Keer, 2000a). 
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3.3 Non-Negative Least Squares (NNLS) 


In this section we show how a QP problem with positive dehnite Hessian 
matrix having the special form (16) can be effectively solved as a nonnegative 
least squares problem. 

Thanks to property (15), matrix H admits a Cholesky factorization H = C'C. 
Hence we can theoretically recast problem (16) as the Non-Negative Least 
Squares (NNLS) problem: 

miup -||Cp — C“^u ||2 (20a) 

s.t. p > 0 (20b) 

A simple and effective active-set method for solving the NNLS problem (20) 
is the one in (Lawson and Hanson, 1974, p.l61), that is extended here in 
Algorithm 3 to directly solve (16) without explicitly computing the Cholesky 
factor C and its inverse and to handle warm starts. After a hnite number 
of steps. Algorithm 3 converges to the optimal contact force vector p* and 
returns the normal displacement vector u* whose components Pij, Uij satisfy 
Pi,j > 0, Uij > Uij, {uij - Uij)pij = 0, and (13), V(i, j) G Ic- 

The method is easy to warm start in case of a loading scenario consisting of an 
alternating sequence of increasing or decreasing far-held displacements. The 
contact forces determined for a given imposed displacement are used to initial¬ 
ize vector p. Due to the monotonicity of the contact solution, this initialization 
is certainly much closer to the optimal solution p* than a zero vector. This 
usually signihcantly reduces the iterations of the method to convergence. Such 
a warm start has a fast implementation requiring a projection of the forces 
of the points belonging to /^(A^) to the same points of the trial domain 
/^(Afc_|_i) for a new imposed far held displacement A^+i. For an increasing 
far-held displacement, i.e., A^+i > A^ the forces in the elements belonging 
to /^(Afc_|_i) — /^(Afc) are simply initialized equal to zero. In the numerical 
experiments of Section 4, the time required for this projection will be added to 
the global solution time for a consistent comparison with the greedy method 
with cold start and with the constrained CG algorithm. 
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Algorithm 3. Non-Negative Least Squares (NNLS) 

Input: Matrix H = H' 0, vector u, initial guess p; maximum number it'max 
of iterations, tolerance e > 0. 

(1) X <(— {i G {1,..., n} : Pi > 0}; init ^ FALSE; /c <(— 0; 

(2) if X = 0 then init ^ TRUE; 

(3) w ^ Hp — u; 

(4) if ((w > — e or X = {1,..., n}) and init = TRUE) or k > Amax then 
go to Step 13; 

(5) if init = TRUE then i ^ argminig{i^,,,^„}\x w*; X <(— X U {i}; 
else init ■(— TRUE; 

(6) s/ •(— solution of the linear system H/S/ = u/ 

(7) if Sx > — e then p ^ s and go to Step 3; 

(8) j ^ argmin;,6X; s,<o 

(9) p^p + ^(s-p); 

(10) Xq {/i G X : p/i = 0}; 

(11) X 4— X \ Xq; fc /c + 1; 

(12) go to Step 6; 

(13) p* ^ p; 

(14) u* ^ w + u; 

(15) end. 

Output; Contact force vector p* and normal displacement vector u* satisfy¬ 
ing u* = Hp, u* > u, p* > 0, (u* - u)'p = 0. 

Note that Step 6 of Algorithm 3 is equivalent to Step 2.4 of Algorithm 1 
and has been performed by using the MATLAB’s mldivide solver. This step 
can be accelerated by the use of an approach based on the EFT (for its im¬ 
plementation, see e.g. (Batrouni et ah, 2002)). Alternatively, since the set Xq 
changes incrementally during the iterations of the algorithm, more efficient it¬ 
erative QR (Lawson and Hanson, 1974, Chap. 24) or LDL^ Bemporad (2014) 
factorization methods can be employed. 


3.3.1 Warm-started NNLS via accelerated Gradient Projection (NNLS-hGP) 

An alternative method to solve Problem (16) is to use an accelerated gradi¬ 
ent projection (GP) method for QP (Nesterov, 1983; Patrinos and Bemporad, 
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Algorithm 4. Accelerated Gradient Projection (GP) 

Input: Matrix H = H' >- 0 and its Frobenius norm L, vector u, initial guess 
p, number K of iterations. 

(1) p ^ p; 

(2) for i = 0,..., JF — 1 do; 

(2.1) /3 = max{^,0}; 

(2.2) s = p + /3(p-p); 

(2.3) w = Hs — u; 

(2.4) p ^ p; 

(2.5) p ^ max{s — ^w, 0}; 

(3) end. 

Output; Warm start for contact force vector p and elastic correction vector 

w. 


2014). Because of the simple nonnegative constraints in (16), rather than go¬ 
ing to the dual QP formulation as in (Patrinos and Bemporad, 2014), we for¬ 
mulate the GP problem directly for the primal QP problem (16). Numerical 
experiments have shown slow convergence of a pure accelerated GP method 
to solve (16). However, we can use the method to warm start Algorithm 3, 
as described in Algorithm 4. If Algorithm 4 is executed {K > 0), it returns a 
vector p that is immediately used as an input to Algorithm 3, otherwise one 
can simply set p = 0 (cold start). As shown in Section 4, GP iterations provide 
large benefits in warm starting the NNLS solver, therefore allowing taking the 
best advantages of the two methods: quickly getting in the neighborhood of 
the optimal solution (GP iterations of Algorithm 4) and getting solutions up 
to machine precision after a hnite number of iterations (the active-set NNLS 
Algorithm 3). 


3-4 Alternating Direction Method of Multipliers (ADMM) 


The QP problem (16) can also be solved by the Alternating Direction Method 
of Multipliers (ADMM), which belongs to the class of augmented Lagrangian 
methods. The reader is referred to (Boyd et ah, 2011) for mathematical details. 
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The method treats the QP (16) as the following problem 


minp,s ip': 


p'Hp - u'p + ^(s) 


( 21 ) 


s.t. p = s 


where 


9{s) 


1 


+00 if s < 0 


0 if s > 0 


Then, the augmented Lagrangian function 



is considered, where p > 0 is a parameter of the algorithm. The basic ADMM 
algorithm consists of the following iterations: 

pfc+i _ argmiup Lp(p, s^, w^) 


( 22 ) 


gfc +1 _ argmiUs Lp(p^’''^, s, w^) 
+ p(p^+^ — s^+^) 


A scaled form with over-relaxation of the ADMM iterations (22) is summarized 
in Algorithm 5. The algorithm is guaranteed to converge asymptotically to the 
solution p*, u* of the problem. The over-relaxation parameter « > 1 is intro¬ 
duced to improve convergence, typical values for a suggested in (Boyd et ah, 
2011) are a e [1.5,1.8]. 

A warm start of the algorithm that takes into account the loading history 
is possible in a way analogous to that described for the NNLS approach of 
Section 3.3. However, as an additional complexity, also an initialization for 
the dual variable vector w must be provided, possibly obtained by projecting 
the solution obtained for a certain to that for A^+i. 
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Algorithm 5. Alternative Direction Method of Multipliers (ADMM) 

Input: Matrix H = H' >- 0, vector u, initial guesses p, w, parameter p > 
0, over-relaxation parameter a > 1, maximum number i^max of iterations, 

tolerance e > 0._ 

(1) M^(iH + I)-^ 

(2) Wp ^ -iw; 

(3) s ^ p; 

(4) z ^ 0; 

(5) while (z < J ^max and ||p — s||oo > e) or z = 0 do; 

(5.1) s ^ M(p - Wp - in); 

(5.2) s ^ os -I- (1 — q;)p; 

(5.3) p ^ max{s-|-Wp, 0}; 

(5.4) Wp ^ Wp -h s - p; 

(5.5) z •(— z -f 1; 

(6) p* ^ p; 

(7) u* ^ u - pwp; 

(8) end. 

Output; Contact force vector p* and normal displacement vector u* satisfy¬ 
ing u* = Hp, u* > u, p* > 0, (u* - u)'p = 0. 

4 Performance comparison of the algorithms 


The optimization algorithms presented in the previous section are herein ap¬ 
plied to the frictionless normal contact problem between a numerically gen¬ 
erated pre-fractal rough surface and a half-plane, in order to compare their 
performance in terms of number of iterations required to achieve convergence 
and computation time. 

The random midpoint displacement algorithm (Peitgen and Saupe, 1988) is 
used to generate the synthetic height field of surfaces with multiscale fractal 
roughness, i.e., with a power spectral density (PSD) function of the height 
field of power-law type. The surface with a given resolution (pre-fractal) is re¬ 
alized by a successive rehnement of an initial coarse representation by adding 
a sequence of intermediate heights whose elevation is extracted from a Gaus¬ 
sian distribution with a suitable rescaled variance, see a qualitative sketch 
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Progressive surface refinement 


Figure 3. Rough surfaces with multi-scale roughness and different resolution, nu¬ 
merically generated by the random midpoint displacement algorithm, 
in Fig. 3. Several applications of the method to model rough surfaces for 
contact mechanics simulations are available in (Zavarise et ah, 2004, 2007; 
Paggi and Ciavarella, 2010). 

In particular, we consider a test problem consisting of a surface with Hurst 
exponent H = 0.7, lateral size L = 100 /im and 512 heights per side, which 
corresponds to the highest discretization used to sample real surfaces with 
a confocal prohlometer, like the Leica DCM3D available at the Multi-scale 
Analysis of Materials (MUSAM) Laboratory of IMT Lucca, Italy. Similar dis¬ 
cretizations are obtained in case of AFM. The surface is brought into contact 
with an elastic half-plane under displacement control. Ten displacement steps 
are imposed to reach a maximum far-held displacement which is set equal to 
(■Cmax — 'Cave)/2, where ^rastx aud .^ave are the maximum and the average ele¬ 
vations of the rough surface, respectively. All the simulations are carried out 
with the server 653745-421 Proliant DL585R07 from Hewlett Packard with 
128 GB Ram, 4 processors AMD Opteron 6282 SE 2.60 GHz with 16 cores 
running MATLAB R2014b. 

The parameters for the Greedy CG method are the maximum number of 
iterations iPmax = 1 x 10® and the convergence tolerance e = 1 x 10“®. The 
contact forces are initialized at zero (cold start). The constrained CG method 
also considers iPmax = 1 x 10® and the same tolerance e = 1 x 10“®. Both the 
original version by Polonsky and Keer (1999) (labeled P&K1999 in Fig. 4) 
and its warm-started variant (labeled P&K1999 + warm start in Fig. 4) are 
considered. 

For the NNLS algorithm (Algorithm 3) we adopt the warm start strategy 
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Dimensionless force, P/(EA) 

Figure 4. Comparison between the optimization algorithms in terms of computation 
time. 

based on the projection of contact forces from the solution corresponding 
to a previous displacement step. Alternatively, for NNLS+GP, 100 gradient 
projections are used to initialize vector p. For the ADMM method we use 
a = 1.5, p = 1, iFmax = 3 X 10^ and e = 10“®. The total number n of 
optimization variables is varying with A and therefore with the force level. 
For the highest indentation we have n = 35555. Warm starting the algorithm 
is achieved by projecting primal variables as for the NNLS and dual variables 
w as well. The projection simply consists of assigning the values of p*j and 
w* of the boundary elements in contact for the step to the same boundary 
elements belonging to the trial contact domain Ic corresponding to the higher 
indentation A^+i. 

Once convergence is achieved for each imposed far-held displacement, the 
optimization algorithms provide the same normal force P and contact do¬ 
mains, with small roundoff errors due to hnite machine precision. The CPU 
time required by each method to achieve convergence are shown in Fig. 4 
vs. the dimensionless normal force P/{EA), where E is the Young’s mod¬ 
ulus and A = is the nominal contact area. The best performance is 
achieved by the application of the NNLS method with 100 gradient projec¬ 
tions (GP), which is 26 times faster than the original constrained CG method 
by Polonsky and Keer (1999) and about two orders of magnitude faster than 
the ADMM and the Greedy CG algorithms. 
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Figure 5. Computation times for the Greedy method for different sizes n of the 
contact superset Ic- CG vs. MATLAB’s mldivide solver and CG vs. Gauss-Seidel 
algorithm. 

As outlined in the introduction, the Greedy method can be used in conjunction 
with other algorithms for solving the unconstrained linear system of equations 
(Step 2.4) than the CG algorithm. Although an extensive comparison of dif¬ 
ferent solvers of linear systems of equations with positive dehnite matrices is 
ontside the scope of this paper, we tested the Greedy algorithm after replacing 
the CG Step 2.4 with the optimized built-in mldivide function of MATLAB, 
or with the Gauss-Seidel algorithm, as proposed in Borri-Brunetto et ah (1999, 
2001 ). 

The MATLAB’s mldivide solver (which employs the Cholesky factorization) 
leads to a reduction of computation time of 30 — 40%, almost regardless of the 
size of the system n, see Fig. 5. Even with this gain in computation speed, 
the overall performance is still quite far from that of the NNLS Algorithm 3 
on the platform used for the tests. Moreover, the MATLAB solver leads to an 
error of lack of memory for n > 20000, a serious problem for large systems 
that is not suffered by the CG solver described in Step 2.4 of Algorithm 1. 
The Gauss-Seidel algorithm does not suffer for the lack of memory but it is 
about 3 times slower than the CG method. 

The effect of the number K of GP iterations applied before the NNLS algo¬ 
rithm is investigated in Fig. 6 for the same test problem whose results were 
shown in Fig. 4. By increasing K from 0 to 100 we observe a reduction in 
the total compntation time due to a decrease in the number of iterations re- 



• mldivide vs. CG 
■ Gauss-Seidel vs. CG 
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P/(EA) 


Figure 6. Computation times of the NNLS algorithm depending on the number K 
of gradient projection (GP) iterations. 

quested by the NNLS algorithm to achieve convergence thanks to a better 
initial guess of p. However, a further increase in K (see, e.g., the blue curve 
in Fig. 6 corresponding to K = 200 iterations) does not correspond to further 
savings of CPU time. This is due to the fact that the number of NNLS itera¬ 
tions was already reduced to its minimum for K = 100 GP iterations, so that 
the application of further gradient projections are just leading to additional 
CPU time without further beneht. 


5 Cascade multi-resolution (CMR) method 


5.1 Algorithm 


A further speed-up of computation time, as compared to the NNLS method, 
can be achieved by improving the criterion for the guess of the initial set Ic of 
points in contact. The standard criterion based on checking the interpenetra¬ 
tion of the surface heights into the half-plane in case of a rigid body motion 
is the most conservative. However, at convergence, only a small subset of 
that initial set is actually in contact. Therefore, a better choice of the initial 
trial contact domain would reduce the size of the system of linear equations 
with an expected beneht in terms of computation time. 
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As shown in (Borri-Brnnetto et al., 1999) via nnmerical simulations on pre¬ 
fractal surfaces with Hurst exponent H > 0.5 and different resolution, by 
rehning the surface height held via a recursive application of the random 
midpoint displacement algorithm the real contact area of each surface repre¬ 
sentation decreases by reducing 5, as illustrated in the sketch in Fig. 7. In the 
fractal limit of 5 —)■ 0, the real contact area vanishes. Therefore, this property 
of lacunarity implies that the heights that are not in contact for a coarser 
surface representation are not expected to come into contact by a successive 
rehning of the height held, for the same imposed far-held displacement. 


Therefore, as a better criterion, the initial trial contact domain can be selected 
by retaining, among all the heights selected by the rigid body interpenetration 
check, only those located within the areas of inhuence of the nodes belonging 
to the contact domain of a coarser representation of the rough surface for the 
same imposed displacement A. 


As graphically shown in Fig. 7, an area of inhuence of a given node in contact 
can be dehned by the radius \/2(5, where 5 is the grid size of the coarser surface 
representation. Since the criterion is not exact, it is convenient to consider a 
multiplicative factor h larger than one for the radius dehning the nodal area of 
inhuence. It is remarkable to note that this numerical scheme can be applied 
recursively to a cascade of coarser representations of the same rough surface. 
As a general trend, computation time is expected to drastically diminish by 
increasing the number of cascade projections. However, the propagation of 
errors due to the wrong exclusion of heights that would actually make contact 
cannot be controlled by the algorithm and it is expected to increase with the 
number of projections as well. The advantage of the method is represented 
by the fact that, in addition to saving computation time with respect to that 
required by the NNLS algorithm to solve just the contact problem for the 
hnest surface, all the contact predictions for the coarser scale representations 
of the same surface will be available for free, which is a useful result for the 
multi-scale characterization of contact problems. Moreover, the CMR method 
can be used in conjunction with any of the optimization algorithms presented 
in the previous sections. 
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Detected by the geometric 
criterion but not in contact 
after the simulation 


-► 

Surface refinement 


Figure 7. A sketch illustrating the property of lacunarity of the contact domain: the 
real contact area progressively diminishes by refining the surface, until vanishes in 
the fractal limit of (5 —)• 0. This implies that some boundary elements detected by 
the rigid-body interpenetration criterion (dashed grey elements) can be neglected a 
priori since they are outside the real contact area corresponding to the coarse scale 
contact solution. 


Algorithm 6. Cascade multi-resolution (CMR) algorithm 

Input; s = 1,... ,l surface representations with different resolution or grid 
spacing 5(s); area of influence parameter h > 1. 

(1) for s = 1,..., I do: 

(1.1) Determine Ic(s) = {(i,j) E In(s) : (ij > ( max(^ )-A}; 

(1.2) if s = 1 then Ic,p(s) = Ic(s) 

else Ic,p(s) = E Ic(s) : = ||xij - < hd(s - 1)}, 

V(k,l)El^(s-l) 

end 

(2) Construct H based on the projected trial contact domain ^c,p(s); 

(3) Apply optimization algorithms (e.g., NNLS) and determine p*, u*, 

(4) end. 


The algorithm is illustrated in Algorithm (6). 

















































5.2 Validation in case of numerically generated and real rough surfaces 


To assess the computational performance of the approach described in Sec¬ 
tion 5.1, the CMD method is applied in conjunction with the NNLS algorithm 
to pre-fractal surfaces with different H numerically generated by the RMD 
method. As an example, the lateral size is 100 pm for all the surfaces and the 
hnest resolution whose contact response has to be sought corresponds to 256 
heights per side. The method requires the storage of the coarser representa¬ 
tions of such surfaces that are in any case available by the RMD algorithm 
during its various steps of random addition. 

We apply the cascade of projections starting by a coarser representation of 
the surfaces with only 16 heights per side and then considering 32, 64, 128 
and hnally 256 heights per side. A parameter h = 2 has been used for the 
dehnition of the area of influence. The solution of the contact problem for 
the surface with 16 heights per side is obtained in an exact form since it is 
the starting point of the cascade, whereas the contact predictions for the hner 
surface representations can be affected by an error intrinsic in the criterion. 
The approximate predictions for the surface with 256 heights per side are 
compared with the reference solution corresponding to the application of the 
NNLS algorithm with warm start directly to the hnest representation of the 
rough surface. 

The computation time of the CMR-I-NNLS solution is the sum of the CPU 
time required to solve all the coarser surface representations and it is found to 
be much less than the CPU time required by the NNLS algorithm to solve just 
one single surface with the hnest resolution, see Fig. 8, where we observe a 
reduction of 50% in CPU time almost regardless of H. The relative error in the 
computation of the maximum normal force between the predicted solution and 
the reference one is a rapid decreasing function of iL, as shown in Fig. 8(d). 
Considering that real surfaces have often a Hurst exponent H > 0.5, this is 
very promising. 

A synthetic diagram illustrating the ehect of the parameter h for the surface 
with H = 0.7 and for a single imposed displacement corresponding to the 
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(a) H = 0.3 (b) H = 0.5 



(c) H = 0.7 (d) Relative error 


Figure 8. Performance of the CMR+NNLS method applied to numerically generated 
fractal surfaces with different Hurst exponent H, h = 2. 

maximum load in Fig. 8(c) is shown in Fig. 9. The relative error is rapidly 
decreasing to values less than 1% by increasing h. The ratio between the 
number of points expected to be in contact after the application of the CMR 
projection criterion, and the number of points that would be included by 
using the classic rigid-body interpenetration check, n, is ranging from 0.4 to 
0.8 by increasing h from 1.25 to 3.0. The ratio between CPU times, on the 
other hand, tends to an asymptotic value of 0.6, which implies a saving of 40% 
of computation time as compared to the exact solution, with less than 0.01% 
of relative error. 

We also check the CMR method for warm starting on real surfaces not dis¬ 
playing the ideal fractal scaling at any length scale, to better assess possible 
limits of applicability. As a practical example we consider the surface of tex¬ 
tured silicon solar cells sampled with two different lenses in order to achieve 
two different magnihcations (lOx and lOOx) by using the confocal prohlometer 
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Figure 9. Performance of the CMR+NNLS method with respect to NNLS for a 
numerically generated fractal surface with H = 0.7, depending on the parameter h. 

Leica DCM3D, see Fig. 10. The PSD function of such a surface sampled with 
512 points per side presents a power-law trend for high frequencies (fine resolu¬ 
tions) and a cut-off to the power-law at low frequencies (coarse resolutions). In 
the power-law regime the surface is characterized by a Hurst exponent H = 0.6 
that can be determined by the slope of the PSD function as customary. 

As a main difference with respect to pre-fractal rough surfaces generated by 
the RMD algorithm, the application of the CMR method requires a hlter to 
downsample the acquired surfaces and extract their coarser representations. 
The CMR method is applied to the two surfaces acquired with lOx and lOOx 
magnihcations using h = 1.5 and considering a cascade of projections in¬ 
volving coarser representations of the finest surfaces with 64 and 128 heights 
per side. A single contact step corresponding to an imposed far-held normal 
displacement equal to (^max — ^ave)/5 is examined. 

The application of the CMR-I-NNLS method to the surface acquired at lOOx 
leads to very good results in line with those observed for ideal fractal surfaces. 
The relative error in the prediction of the normal load is —0.4%, with a sav¬ 
ing of CPU time of 18% as compared to the direct application of the NNLS 
algorithm. On the other hand, the method applied to the surface acquired at 
lOx leads to poor results in terms of accuracy with —98% of relative error and 
almost no saving in computation time. This bad performance is due to the fact 
that the property of lacunarity of the contact domain, strictly connected with 
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0 0 0 0 


(a) lOx (b) lOOx 

Figure 10. Surface of textured silicon solar cells sampled with a confocal profilometer 
at two different magnifications (lOx and lOOx) obtained by using two different lenses. 



Figure 11. Power spectral density function (PSD) of the two sampled rough surfaces 
shown in Fig. 10. 


the self-affine scaling of roughness due to fractality, does not hold anymore 
for the surface sampled at lOx due to the cut-off to its power-law PSD. As 
a consequence, the CMR method erroneously excludes many possible points 
from the initial contact domain suggested by the rigid body interpenetration 
check that are actually relevant for contact. Therefore, in conclusion, the CMR 
method is efficient for warm starting the NNLS algorithm, but it should be 
strictly applied to numerically generated or real rough surfaces provided that 
the self-affine properties of roughness are conhrmed by a PSD function of 
power-law type. 
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6 Conclusion 


This paper has shown how the problem of frictionless normal contact between 
rough surfaces within the BEM framework can be solved very efficiently by 
exploiting ideas from convex quadratic programming. A series of efficient op¬ 
timization algorithms has been proposed and compared with the traditional 
Greedy method and constrained CG algorithm. As the lack of convergence 
of the Greedy method seems to be a rare phenomenon, it remains an open 
question to establish the conditions on H and u for which the algorithm is 
guaranteed to converge. 

The NNLS algorithm warm started by accelerated gradient projections was 
shown at least two orders of magnitude faster than the Greedy method and 
26 times faster than the original constrained GG algorithm. 

Finally, we explored another method for warm starting the optimization al¬ 
gorithms, this time focusing on a selective reduction of the size of the initial 
trial contact domain based on the multi-resolution properties of roughness. 
The resulting cascade multi-resolution (GMR) method allows a further sav¬ 
ing of about 50% of GPU time as compared to NNLS for contact simulations 
involving numerically generated fractal surfaces. Relative errors were found 
less than 2% for surfaces with H > 0.5, by using h = 2, that was found a 
good compromise between accuracy and computation time. Moreover, it has 
to be remarked that not only the solution of the finest contact problem is 
gained by the GMR-I-NNLS method with much less GPU time, but also the 
contact problems involving all the coarser representations of the finest sur¬ 
face. These results are particularly important for speeding up intensive Monte 
Garlo simulations involving a sequence of contact simulations for a population 
of fractal surface with different resolution. So far, to the best of the authors’ 
knowledge, such extensive simulations, that are important to determine more 
reliable trends from the statistical point of view, have been limited to popu¬ 
lations of 20 to 50 randomly generated surfaces. 

In case of real surfaces, a very good performance (less than 2% of error with 
3 cascades and at least 18% of GPU time saved for one single imposed dis- 
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placement step) has been demonstrated in case of power-law PSDs, assuring 
the self-affine scaling of roughness which represents the main underlying as¬ 
sumption for the algorithm applicability. For surfaces with a cut-off to the 
power-law PSD, on the other hand, the CMR-I-NNLS method has given poor 
results in terms of accuracy and in any case almost no saving in CPU time as 
compared to the pure application of NNLS. Therefore, this warm start method 
should be used with care and only in a range where the PSD is of power-law 
type. 

Finally, we point out that the proposed optimization methods can also be 
applied to frictional contact problems by using for instance the complete BEM 
formulation as in (Pohrt and Li, 2014). Although this issue is left for further 
investigation, we expect an even more signihcant gain in CPU time by applying 
the algorithms presented in this paper instead of other optimization methods, 
since the size of the problem is by far signihcantly increased as compared to 
the frictionless case. 
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